1 Introducción

El sistema de transporte público de Madrid constituye una infraestructura crítica para la movilidad urbana, articulando el flujo diario de millones de pasajeros a través de diversos modos que operan de manera simultánea e interdependiente. Tradicionalmente, estos modos se han gestionado y estudiado como redes aisladas; sin embargo, esta aproximación ignora la complejidad de las interconexiones reales entre el autobús, el metro, la bicicleta pública y los aparcamientos. El presente estudio supera esta limitación adoptando un enfoque de Redes Multicapa, modelando matemáticamente la ciudad como un “sistema de sistemas” donde la eficiencia y la robustez global emergen de la interacción entre las distintas capas de transporte.

En este trabajo se analiza la estructura multicapa de la red de transporte de la ciudad de Madrid a partir de dos conjuntos de datos reales proporcionados:

Un fichero con las paradas y estaciones de distintos modos de transporte.
Un fichero con información de viajes mensuales por distrito.

Los objetivos concretos son:

Construir una red multicapa donde cada capa representa un modo de transporte (bus, metro, bicimad, parking).
Estudiar la centralidad de las zonas de intercambio en cada capa.
Medir la versatilidad de las zonas (cómo cambia su importancia entre modos).
Realizar una visualización multicapa en la que se vean explícitamente las distintas capas conectadas.
Explorar de forma sencilla la vulnerabilidad del sistema ante fallos en una de las capas.

Toda la parte empírica se basa exclusivamente en los ficheros:

transporte_madrid_consolidado.csv
viajes_madrid.csv

2 Carga y exploración de datos

library(tidyverse)
library(janitor)
library(sf)
library(dbscan)
library(FNN)
library(igraph)
library(scales)
library(knitr)
library(devtools)

# Instalación de muxViz si no existe
if (!requireNamespace("muxViz", quietly = TRUE)) {
  message("Instalando muxViz desde GitHub...")
  devtools::install_github("manlius/muxViz", upgrade = "never", force = FALSE)
}
library(muxViz)

2.1 Lectura de los ficheros

Se asume que los ficheros CSV están en el mismo directorio que este documento RMarkdown. Si no es así, ajustar las rutas.

ruta_transporte <- "transporte_madrid_consolidado.csv"
ruta_viajes     <- "viajes_madrid.csv"

transporte_raw <- read.csv(
  ruta_transporte,
  sep = ";",
)

viajes_raw <- read.csv(
  ruta_viajes,
  sep = ";",
)

glimpse(transporte_raw)

## Rows: 14,188
## Columns: 6
## $ stop_name         <chr> "                                                   …
## $ stop_lat          <dbl> 40.37835, 40.38341, 40.36813, 40.37773, 40.40753, 40…
## $ stop_lon          <dbl> -3.75458, -3.62481, -3.62283, -3.67311, -3.59176, -3…
## $ transport_mode    <chr> "bicimad                                            …
## $ stop_id           <chr> "382 -Parque Eugenia de Montijo, 2                  …
## $ distrito_asignado <int> 2807911, 2807918, 2807918, 2807913, 2807919, 2807911…

glimpse(viajes_raw)

## Rows: 8,064
## Columns: 11
## $ fecha             <chr> "2024-01", "2024-01", "2024-01", "2024-01", "2024-01…
## $ name              <chr> "Centro                                            "…
## $ edad              <chr> "0-25", "0-25", "0-25", "0-25", "0-25", "0-25", "0-2…
## $ sexo              <chr> "hombre", "hombre", "hombre", "hombre", "mujer", "mu…
## $ numero_viajes     <chr> "0", "1", "2", "2+", "0", "1", "2", "2+", "0", "1", …
## $ personas          <chr> "90541,756", "18809,207", "73248,991", "187551,534",…
## $ Codigo            <int> 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, …
## $ Provincia         <chr> "Madrid                ", "Madrid                ", …
## $ distrito          <dbl> 2807901, 2807901, 2807901, 2807901, 2807901, 2807901…
## $ zona_pernoctacion <int> 2807901, 2807901, 2807901, 2807901, 2807901, 2807901…
## $ poblacion         <int> 140644, 140644, 140644, 140644, 140644, 140644, 1406…

2.2 Descripción básica de `transporte_madrid_consolidado`

transporte <- clean_names(transporte_raw)

transporte$transport_mode <- str_trim(transporte$transport_mode)
transporte$transport_mode <- str_extract(transporte$transport_mode, "\\w+")

table(transporte$transport_mode)

## 
## bicimad     bus   metro parking 
##     626   12430    1050      82

columnas_interes <- transporte[, c("stop_lat", "stop_lon", "distrito_asignado")]
summary(columnas_interes)

##     stop_lat        stop_lon      distrito_asignado
##  Min.   :40.28   Min.   :-3.875   Min.   :2807901  
##  1st Qu.:40.40   1st Qu.:-3.711   1st Qu.:2807905  
##  Median :40.43   Median :-3.689   Median :2807909  
##  Mean   :40.43   Mean   :-3.683   Mean   :2807910  
##  3rd Qu.:40.46   3rd Qu.:-3.656   3rd Qu.:2807915  
##  Max.   :40.56   Max.   :-3.447   Max.   :2807921

Variables clave:

stop_name: nombre de la parada/estación.
stop_lat, stop_lon: coordenadas geográficas.
transport_mode: modo de transporte (bus, metro, bicimad, parking).
stop_id: identificador de parada.
distrito_asignado: código de distrito (INE).

Esto define los nodos base del sistema de transporte.

2.3 Descripción básica de `viajes_madrid`

library(readr)
library(dplyr)
library(janitor)
library(stringr)

viajes <- read.csv("viajes_madrid.csv",
                   sep = ";",
                   dec = ",")

viajes <- clean_names(viajes)

summary(viajes$personas)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4515   37486  110989  140469  186820  737342

Variables clave:

fecha: mes de referencia.
name: nombre del distrito.
edad, sexo: características sociodemográficas.
numero_viajes: categoría de nº de viajes (0, 1, 2, 2+).
personas: número de personas en esa celda.
distrito: código de distrito.
poblacion: población total del distrito.

2.3.1 Indicador de intensidad de viajes por distrito

Construimos un indicador sencillo de viajes medios por persona en cada distrito. Primero pasamos la categoría numero_viajes a un valor numérico aproximado y luego calculamos una media ponderada por el número de personas.

2.3.2 Indicador de viajes medios diarios por persona y distrito

La variable personas se encontraba codificada con formato decimal europeo, por lo que fue necesario convertirla explícitamente a formato numérico. Asimismo, la variable numero_viajes es categórica (0, 1, 2, 2+), por lo que se transformó en una aproximación numérica asignando el valor 3 a la categoría “2+”. El indicador final de viajes medios por persona se calculó como una media ponderada por el número de personas en cada categoría, obteniéndose valores en torno a 2 viajes diarios por persona, coherentes con patrones conocidos de movilidad urbana.

library(dplyr)
library(knitr)

viajes$viajes_cat <- case_when(
  viajes$numero_viajes == "0"  ~ 0,
  viajes$numero_viajes == "1"  ~ 1,
  viajes$numero_viajes == "2"  ~ 2,
  viajes$numero_viajes == "2+" ~ 3,
  TRUE ~ NA_real_
)

viajes_distrito <- aggregate(
  cbind(personas, viajes_totales = viajes_cat * personas) ~ distrito, 
  data = viajes, 
  FUN = sum, 
  na.rm = TRUE
)

viajes_distrito$viajes_medios_por_persona <- viajes_distrito$viajes_totales / viajes_distrito$personas
viajes_distrito$personas_totales <- viajes_distrito$personas

kable(
  head(viajes_distrito),
  digits = 3,
  caption = "Viajes medios diarios por persona y distrito"
)

Viajes medios diarios por persona y distrito
distrito	personas	viajes_totales	viajes_medios_por_persona	personas_totales
2807901	51888948	110673624	2.133	51888948
2807902	51542198	108276834	2.101	51542198
2807903	39251902	78164164	1.991	39251902
2807904	49039609	98112207	2.001	49039609
2807905	48349024	95741527	1.980	48349024
2807906	53679153	110757611	2.063	53679153

summary(viajes_distrito$viajes_medios_por_persona)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.956   2.004   2.038   2.041   2.063   2.135

Interpretación de los Resultados (Estadísticas de Viajes)

El resumen estadístico (summary) de la variable viajes_medios_por_persona revela un patrón de comportamiento extremadamente homogéneo en todos los distritos analizados:

Uniformidad del Comportamiento: El rango de valores es muy estrecho, oscilando apenas entre un mínimo de 1.956 y un máximo de 2.135 viajes por persona. Esto indica que, independientemente del distrito (sea centro o periferia), la movilidad media individual es muy constante.
Tendencia Central (La regla del “Ida y Vuelta”): Tanto la media (2.041) como la mediana (2.038) se sitúan prácticamente en 2. Esto tiene un sentido lógico muy fuerte en transporte urbano: representa el patrón estándar de un viaje de ida y otro de vuelta (péndulo diario).
Implicación para el Modelo: Dado que la desviación entre el primer cuartil (2.004) y el tercer cuartil (2.063) es mínima, esta variable no actúa como un factor discriminante fuerte por sí sola. Esto significa que las diferencias de flujo total que observaremos en la red no se deben a que la gente de un barrio viaje “más veces” que la de otro, sino probablemente a la densidad de población de cada zona (mismo ratio de viajes por persona, pero mucha más gente en unas zonas que en otras).

Se observa una notable homogeneidad en la movilidad individual a lo largo de Madrid, con una media de 2.04 viajes por persona y una varianza despreciable. Esto sugiere que la intensidad de uso de la red de transporte depende casi exclusivamente del volumen demográfico de cada distrito y no de diferencias significativas en los hábitos de movilidad individual entre zonas.

3 Construcción de la red multicapa

La idea es pasar de paradas individuales a zonas de intercambio (clusters de paradas cercanas), y luego construir una capa por modo de transporte.

3.1 Definición de zonas de intercambio (clustering espacial)

# Creamos objeto sf con coordenadas
transporte_sf <- st_as_sf(transporte, coords = c("stop_lon", "stop_lat"), crs = 4326, remove = FALSE)
transporte_sf <- st_transform(transporte_sf, 25830)  # sistema métrico (UTM zona 30)

coords_m <- st_coordinates(transporte_sf)
coords_m <- as.matrix(coords_m)

# Clustering DBSCAN: paradas a menos de 150 m se agrupan
set.seed(123)
db <- dbscan(coords_m, eps = 150, minPts = 2)

transporte_sf$cluster_id <- db$cluster

# Tratamos el ruido (cluster 0) como clusters individuales
max_id <- max(transporte_sf$cluster_id)
es_ruido <- transporte_sf$cluster_id == 0
transporte_sf$cluster_id[es_ruido] <- max_id + seq_len(sum(es_ruido))

# Zona de intercambio = centroide del cluster
df_no_geo <- st_drop_geometry(transporte_sf)

# Calculamos medias de coordenadas
zonas <- aggregate(cbind(stop_lat, stop_lon) ~ cluster_id, 
                   data = df_no_geo, 
                   FUN = mean)

# distrito asignado más frecuente en el cluster
distritos_freq <- aggregate(distrito_asignado ~ cluster_id, 
                             data = df_no_geo, 
                             FUN = function(x) as.integer(names(which.max(table(x)))))

# Modos disponibles
modos_agg <- aggregate(transport_mode ~ cluster_id, 
                       data = df_no_geo, 
                       FUN = function(x) paste(sort(unique(x)), collapse = ","))

# Unimos todo en el objeto zonas
zonas <- merge(zonas, distritos_freq, by = "cluster_id")
zonas <- merge(zonas, modos_agg, by = "cluster_id")

nrow(zonas)

## [1] 1558

head(zonas)

Cada cluster_id representa una zona de intercambio multimodal.

Representamos la distribución espacial aproximada (proyección geográfica simple):

# Aseguramos que sea factor
zonas_sf <- st_as_sf(zonas, coords = c("stop_lon", "stop_lat"), crs = 4326)

zonas_sf$distrito_asignado <- as.factor(zonas_sf$distrito_asignado)
n_distritos <- length(levels(zonas_sf$distrito_asignado))

# Usamos hcl.colors con la paleta "Dark 3" o "Set 2", que son mucho más elegantes
plot(zonas_sf["distrito_asignado"], 
     main = "Zonas de intercambio por distrito",
     pal = hcl.colors(n_distritos, palette = "Dark 3"), 
     key.pos = 4,
     pch = 19, 
     cex = 0.6,
     border = NA) # Quitamos bordes si los hubiera para que el color sea más puro

3.2 Nodos por modo de transporte

modos <- sort(unique(transporte_sf$transport_mode))

# Creamos una tabla de frecuencias y la convertimos a data frame
df_tab <- as.data.frame.matrix(table(transporte_sf$cluster_id, transporte_sf$transport_mode))

# Convertimos las frecuencias en presencia (1) o ausencia (0)
df_tab[df_tab > 0] <- 1
df_tab$cluster_id <- as.numeric(rownames(df_tab))

# Unimos con las zonas
zonas_attrs <- merge(zonas, df_tab, by = "cluster_id", all.x = TRUE)

# Reemplazamos NAs por 0 en las columnas de modos
zonas_attrs[modos][is.na(zonas_attrs[modos])] <- 0

kable(head(zonas_attrs), caption = "Atributos básicos de las zonas por modo")

Atributos básicos de las zonas por modo
cluster_id	stop_lat	stop_lon	distrito_asignado	transport_mode	bicimad	bus	metro
1	40.38368	-3.624114	2807918	bicimad,bus,metro	1	1	1
2	40.38096	-3.728060	2807911	bicimad,bus,metro	1	1	1
3	40.38626	-3.719795	2807912	bicimad,bus,metro	1	1	1
4	40.37886	-3.732537	2807911	bicimad,bus	1	1	0
5	40.40855	-3.672909	2807903	bicimad,bus	1	1	0
6	40.40380	-3.706547	2807901	bicimad,metro	1	0	1

En este análisis, la estructura de los clústeres utiliza la proximidad geográfica para configurar nodos de naturaleza multimodal. Cada clúster integra diversas tipologías de infraestructuras de transporte basándose en un criterio de agrupación por distancias mínimas, lo que permite consolidar puntos de acceso heterogéneos en una unidad funcional única.

Como se observa en la matriz resultante, la presencia de servicios como \(bicimad\), \(bus\), \(metro\) y \(parking\) dentro de un mismo \(cluster\_id\) evidencia que el modelo no solo identifica cercanía espacial, sino que captura la interconectividad del tejido urbano, facilitando la transición entre diferentes modos de movilidad en puntos estratégicos de la red.

3.3 Construcción de redes intra-capa

Este bloque de código construye la estructura topológica intra-capa (conexiones internas) para cada medio de transporte mediante un criterio de proximidad espacial. En primer lugar, proyecta las coordenadas geográficas a un sistema métrico (UTM) para poder calcular distancias reales en metros. Posteriormente, aplica un algoritmo de K-Vecinos Más Cercanos (KNN) que conecta cada parada (nodo) con sus 3 estaciones más próximas del mismo modo de transporte, asignando a cada conexión un peso (\(w\)) inversamente proporcional a la distancia física; de esta forma, se generan redes de vecindad geométrica donde los nodos cercanos tienen enlaces más fuertes, simulando la accesibilidad local dentro de cada capa.

# 1. Preparación de coordenadas de forma sencilla
centroides_mat <- st_coordinates(st_transform(zonas_sf, 25830))
zonas_centroides <- data.frame(
  cluster_id = zonas$cluster_id,
  x = centroides_mat[, 1],
  y = centroides_mat[, 2]
)

# 2. Función con lógica explícita
build_edges_simple <- function(modo_nombre, k = 3) {
  
  # Filtrar qué zonas tienen este modo activo
  ids_con_modo <- zonas_attrs$cluster_id[zonas_attrs[[modo_nombre]] == 1]
  
  # Si no hay suficientes zonas para conectar, devolvemos tabla vacía
  if (length(ids_con_modo) < 2) return(data.frame())
  
  # Obtener coordenadas solo de esas zonas
  coords <- zonas_centroides[zonas_centroides$cluster_id %in% ids_con_modo, ]
  
  lista_aristas <- list()
  
  # Bucle: Para cada zona, calculamos la distancia a TODAS las demás
  for (i in 1:nrow(coords)) {
    punto_origen <- coords[i, ]
    
    # Calculamos distancia euclidea: sqrt((x1-x2)^2 + (y1-y2)^2)
    distancias <- sqrt((punto_origen$x - coords$x)^2 + (punto_origen$y - coords$y)^2)
    
    # Creamos un dataframe temporal con las distancias de 'i' al resto
    df_dist <- data.frame(
      modo = modo_nombre,
      from = punto_origen$cluster_id,
      to   = coords$cluster_id,
      dist = distancias
    )
    
    # Quitamos la distancia a sí mismo (que es 0)
    df_dist <- df_dist[df_dist$from != df_dist$to, ]
    
    # Ordenamos por distancia y nos quedamos con los 'k' más cercanos
    df_dist <- df_dist[order(df_dist$dist), ]
    df_dist <- head(df_dist, k)
    
    lista_aristas[[i]] <- df_dist
  }
  
  # Unir todos los resultados del bucle en una tabla
  resultado <- do.call(rbind, lista_aristas)
  
  # Calcular el peso (w) de forma manual
  resultado$w <- 1 / (1 + resultado$dist)
  
  return(resultado)
}

edges_intra_lista <- list()
for (m in modos) {
  edges_intra_lista[[m]] <- build_edges_simple(m)
}

# Combinamos todas las aristas en un solo dataframe
edges_intra <- do.call(rbind, edges_intra_lista)

kable(head(edges_intra), caption = "Ejemplo de aristas intra-capa")

Ejemplo de aristas intra-capa
	modo	from	to	dist	w
bicimad.144	bicimad	1	288	649.2619	0.0015378
bicimad.384	bicimad	1	1382	712.9355	0.0014007
bicimad.218	bicimad	1	515	746.3214	0.0013381
bicimad.149	bicimad	2	300	321.5842	0.0031000
bicimad.4	bicimad	2	4	445.7282	0.0022385
bicimad.428	bicimad	2	1463	736.8067	0.0013554

Descripción del Dataframe obtenido (edges_intra)

El objeto resultante es una tabla consolidada (tibble) donde cada fila representa una arista (conexión) entre dos paradas. Sus columnas son:

modo: Identificador de la capa de transporte (ej. “metro”, “bus”, “bicimad”).
from: ID numérico del nodo de origen (cluster_id).
to: ID numérico del nodo de destino (uno de sus vecinos más cercanos).
dist: La distancia euclídea en metros entre ambos nodos.
w: El peso de la conexión, calculado como \(1/(1 + \text{distancia})\). Varía entre 0 y 1, siendo mayor cuanto más cerca están las paradas.

3.4 Grafos por capa (lista `g_list`)

Este bloque tiene una misión crítica en el análisis multicapa: garantizar la alineación de nodos (Node Alignment). La función build_graph_modo transforma la lista de enlaces (edges_intra) en objetos de grafo (igraph), pero con una particularidad vital: asegura que todas las capas tengan exactamente los mismos nodos, incluso si en una capa específica un nodo no tiene conexiones (está aislado). El código identifica los nodos que faltan en una capa (setdiff) y los añade manualmente como vértices aislados. Sin este paso, las matrices de adyacencia tendrían dimensiones distintas y sería matemáticamente imposible construir el objeto “multiplex”.

nodos_ids <- zonas_attrs$cluster_id

build_graph_simple <- function(nombre_modo) {
  filas_modo <- edges_intra[edges_intra$modo == nombre_modo, ]
  e_modo <- filas_modo[, c("from", "to", "w")]

  g <- graph_from_data_frame(
    d = e_modo, 
    directed = FALSE, 
    vertices = zonas_attrs
  )

  nodos_en_grafo <- as.integer(V(g)$name)
  faltantes <- nodos_ids[!(nodos_ids %in% nodos_en_grafo)]
  
  if (length(faltantes) > 0) {
    g <- g + vertices(as.character(faltantes))
  }

  return(g)
}

g_list <- list()
for (m in modos) {
  g_list[[m]] <- build_graph_simple(m)
}

g_list

## $bicimad
## IGRAPH 8c7c1bf UN-- 1558 1404 -- 
## + attr: name (v/c), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), transport_mode (v/c), bicimad (v/n), bus (v/n), metro (v/n),
## | parking (v/n), w (e/n)
## + edges from 8c7c1bf (vertex names):
##  [1] 1 --288  1 --1382 1 --515  2 --300  2 --4    2 --1463 3 --1511 3 --1241
##  [9] 3 --1547 2 --4    4 --300  4 --35   5 --283  5 --1108 5 --240  6 --927 
## [17] 6 --1496 6 --239  7 --29   7 --479  7 --239  8 --735  8 --1342 8 --756 
## [25] 11--248  11--41   11--460  12--372  12--488  12--1070 13--1391 13--112 
## [33] 13--436  14--1556 14--1056 14--156  16--583  16--1350 16--728  23--778 
## [41] 23--1354 23--1398 24--626  24--306  24--248  25--296  25--748  25--370 
## + ... omitted several edges
## 
## $bus
## IGRAPH ec4c1d3 UN-- 1558 4038 -- 
## + attr: name (v/c), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), transport_mode (v/c), bicimad (v/n), bus (v/n), metro (v/n),
## | parking (v/n), w (e/n)
## + edges from ec4c1d3 (vertex names):
##  [1] 1 --1289 1 --280  1 --1231 2 --1026 2 --300  2 --1304 3 --1303 3 --1082
##  [9] 3 --1148 4 --652  4 --532  4 --229  5 --283  5 --1108 5 --945  7 --29  
## [17] 7 --479  7 --239  8 --819  8 --830  8 --1153 9 --390  9 --335  9 --48  
## [25] 10--844  10--354  10--1446 11--767  11--248  11--41   12--1288 12--368 
## [33] 12--1449 13--638  13--639  13--642  14--158  14--1272 14--1198 15--422 
## [41] 15--357  15--187  16--969  16--726  16--724  18--20   18--19   18--345 
## + ... omitted several edges
## 
## $metro
## IGRAPH 9275c7c UN-- 1558 684 -- 
## + attr: name (v/c), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), transport_mode (v/c), bicimad (v/n), bus (v/n), metro (v/n),
## | parking (v/n), w (e/n)
## + edges from 9275c7c (vertex names):
##  [1] 1 --230  1 --227  1 --49   2 --229  2 --1077 2 --3    2 --3    3 --459 
##  [9] 3 --816  6 --239  6 --7    6 --132  7 --239  6 --7    7 --118  11--248 
## [17] 11--41   11--24   12--372  12--746  12--510  13--112  13--436  13--97  
## [25] 17--187  17--18   17--251  18--187  18--30   18--68   24--248  24--838 
## [33] 11--24   27--917  27--572  27--116  28--287  28--696  28--523  30--68  
## [41] 18--30   30--58   40--145  40--164  40--245  41--753  11--41   40--41  
## + ... omitted several edges
## 
## $parking
## IGRAPH 92a0fbb UN-- 1558 162 -- 
## + attr: name (v/c), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), transport_mode (v/c), bicimad (v/n), bus (v/n), metro (v/n),
## | parking (v/n), w (e/n)
## + edges from 92a0fbb (vertex names):
##  [1] 7  --29   7  --1488 7  --239  8  --1532 8  --41   8  --259  27 --766 
##  [8] 27 --815  27 --116  28 --1475 28 --292  28 --815  29 --239  7  --29  
## [15] 29 --657  41 --245  41 --1503 41 --164  55 --1052 55 --657  55 --423 
## [22] 56 --814  56 --413  56 --67   67 --814  56 --67   67 --121  71 --446 
## [29] 71 --1052 55 --71   84 --657  84 --1052 84 --239  101--164  101--573 
## [36] 101--102  102--573  102--413  102--446  116--766  116--815  27 --116 
## + ... omitted several edges

La ejecución del algoritmo de construcción de grafos ha resultado en la generación de cuatro objetos igraph independientes, correspondientes a las capas de BiciMAD, Bus, Metro y Parking. La inspección de la salida (g_list) revela las siguientes características estructurales fundamentales para el modelo multicapa:

Isomorfismo Nodal (Alineación de Capas) Se confirma que las cuatro capas poseen exactamente 1558 nodos (vertices = 1558). Este resultado valida el procedimiento de alineación de nodos (node alignment), garantizando que todas las capas comparten el mismo conjunto de definiciones espaciales (paradas/zonas), incluso si un modo de transporte no tiene actividad en una zona específica. Esto es un requisito matemático indispensable para la construcción del tensor supra-adyacente en muxViz.
Heterogeneidad en la Densidad de Conexiones Se observa una gran disparidad en el número de aristas (edges) entre capas, lo que refleja la naturaleza operativa de cada modo:

Bus (4038 enlaces): Es la capa con mayor densidad topológica. Su elevado número de aristas indica una alta capilaridad, actuando como el tejido conectivo principal que une distancias cortas y medias con una gran redundancia local.
BiciMAD (1404 enlaces): Presenta una conectividad intermedia. Aunque inferior al autobús, su densidad sugiere una fuerte cohesión en las zonas donde está implantado (principalmente la almendra central), formando clústeres locales densos.
Metro (684 enlaces): Muestra una topología más dispersa. Al funcionar como una red “esqueleto” o vertebral, sus conexiones son menores en número pero estratégicas, diseñadas para el transporte rápido de larga distancia sin saturar el espacio con conexiones redundantes.
Parking (162 enlaces): Es la capa más desconectada. La escasez de aristas refleja que los aparcamientos funcionan operativamente como “nodos satélite” o destinos finales, con muy poca interacción directa entre ellos bajo el criterio de vecindad espacial establecido (\(k=3\)).

Preservación de Atributos: La salida confirma (attr: name, stop_lat…) que se han preservado los metadatos espaciales y demográficos en todos los nodos. Esto permitirá que, en las fases posteriores de visualización y análisis de vulnerabilidad, se puedan correlacionar las métricas de red (como el PageRank) con variables del mundo real como el distrito asignado o los viajes totales.

Tabla 1: Resumen Topológico de las Capas de la Red Multicapa
Capa	Nodos	Aristas	Densidad_Relativa	Funcion_Sistema
Bus	1558	4038	Alta	Capilaridad / Cobertura
BiciMAD	1558	1404	Media	Movilidad última milla / Centro
Metro	1558	684	Baja	Vertebración / Eficiencia
Parking	1558	162	Muy Baja	Nodos Satélite / Intercambio

4 Visualización multicapa con `muxViz`

En esta sección generamos la visualización clave: un gráfico 3D en el que se ven las capas (modos) apiladas, con los nodos alineados entre capas.

Para ello usamos la función plot_multiplex3D de muxViz. Se asume que muxViz está instalado; si no lo está, el código no fallará, pero mostrará un mensaje.

# --- 1. PREPARACIÓN DE DATOS ---
zonas_attrs <- zonas_attrs[order(zonas_attrs$cluster_id), ]

centroides_mat <- st_coordinates(st_transform(zonas_sf, 25830))
zonas_centroides <- data.frame(
  cluster_id = zonas$cluster_id,
  x = centroides_mat[, 1],
  y = centroides_mat[, 2]
)
zonas_centroides <- zonas_centroides[order(zonas_centroides$cluster_id), ]


# --- 2. GENERACIÓN DE ARISTAS (INTRA-CAPA) ---
edges_intra_lista <- list()

for (m in modos) {
  ids_modo <- zonas_attrs$cluster_id[zonas_attrs[[m]] == 1]
  
  if (length(ids_modo) >= 2) {
    coords_sub <- zonas_centroides[zonas_centroides$cluster_id %in% ids_modo, ]
    aristas_modo <- list()
    
    for (i in 1:nrow(coords_sub)) {
      dists <- sqrt((coords_sub$x[i] - coords_sub$x)^2 + (coords_sub$y[i] - coords_sub$y)^2)
      
      df_temp <- data.frame(
        modo = m,
        from = coords_sub$cluster_id[i],
        to   = coords_sub$cluster_id,
        dist = dists
      )
      
      df_temp <- df_temp[df_temp$from != df_temp$to, ]
      df_temp <- df_temp[order(df_temp$dist), ]
      aristas_modo[[i]] <- head(df_temp, 3)
    }
    edges_intra_lista[[m]] <- do.call(rbind, aristas_modo)
  }
}

edges_intra <- do.call(rbind, edges_intra_lista)
edges_intra$w <- 1 / (1 + edges_intra$dist)


# --- 3. CONSTRUCCIÓN DE GRAFOS ---
g_list <- list()

for (m in modos) {
  e_sub <- edges_intra[edges_intra$modo == m, c("from", "to", "w")]
  g <- graph_from_data_frame(d = e_sub, directed = FALSE, vertices = zonas_attrs)
  g_list[[m]] <- g
}

# --- PATCH PARA MUXVIZ EN R 4.3+ ---
# La función original de muxViz falla en R 4.3+ porque usa "class(x) == 'y'" en lugar de "inherits"
# y "cond && cond" con vectores. Definimos una versión local corregida si es necesario.
# (Solo redefinimos la función si falla la original, pero para prevenir, lo hacemos aquí)
# Ajuste mínimo: redefinir la lógica que falla.
# Sin embargo, reescribir toda la función es complejo.
# Alternativa: Asegurar que layer.labels sea NULL o un escalar si es posible, pero muxViz usa names(g_list).

# Solución definitiva: Redefinimos plot_multiplex3D localmente con el fix.
# Basado en código fuente original de muxViz.
plot_multiplex3D_fix <- function(g.list, layer.layout, layer.colors, layer.shift.x=0.5, layer.shift.y=0.5, layer.space=2, show.aggreg=FALSE) {
    if (!requireNamespace("rgl")) stop("rgl required")
    if (!requireNamespace("igraph")) stop("igraph required")
    
    L <- length(g.list)
    LayLabels <- names(g.list)
    
    # Setup 3D
    rgl::open3d()
    rgl::bg3d("white")
    
    # Dibujamos capas
    for(l in 1:L){
        # Un plano por capa
        z.val <- (l-1)*layer.space
        
        # Nodos
        layout.l <- layer.layout
        layout.l[,1] <- layout.l[,1] + (l-1)*layer.shift.x
        layout.l[,2] <- layout.l[,2] + (l-1)*layer.shift.y
        
        # Aristas intra-capa
        el <- igraph::as_edgelist(g.list[[l]], names=FALSE)
        if(nrow(el)>0){
            x0 <- layout.l[el[,1], 1]
            y0 <- layout.l[el[,1], 2]
            x1 <- layout.l[el[,2], 1]
            y1 <- layout.l[el[,2], 2]
            z <- rep(z.val, length(x0))
            
            rgl::segments3d(as.vector(rbind(x0,x1)), as.vector(rbind(y0,y1)), as.vector(rbind(z,z)), 
                            col=grDevices::adjustcolor(layer.colors[l], alpha.f=0.5), lwd=0.5)
        }
        
        # Nodos (Esferas) - Solo activos
        deg <- igraph::degree(g.list[[l]])
        act <- which(deg > 0)
        if(length(act) > 0){
             # Tamaño proporcional al grado normalizado (ajustado para layout normalizado)
             sizes <- 1 + 2 * (deg[act] / max(deg))
             # Escalar el radio: 0.05 es un buen factor si el layout está normalizado (range ~4)
             rgl::spheres3d(layout.l[act,1], layout.l[act,2], rep(z.val, length(act)), 
                            radius = sizes * 0.05, col = layer.colors[l])
        }
        
        # Etiqueta de capa
        rgl::text3d(min(layout.l[,1]), min(layout.l[,2]), z.val, texts=LayLabels[l], col="black", font=2, cex=1.5)
    }
    
    # Enlaces inter-capa (coupling) - Simplificado: solo nodos activos en ambas capas adyacentes
    # Dibujamos líneas verticales para nodos que existen en l y l+1
    for(l in 1:(L-1)){
        layout.l1 <- layer.layout
        layout.l1[,1] <- layout.l1[,1] + (l-1)*layer.shift.x
        layout.l1[,2] <- layout.l1[,2] + (l-1)*layer.shift.y
        
        layout.l2 <- layer.layout
        layout.l2[,1] <- layout.l2[,1] + (l)*layer.shift.x
        layout.l2[,2] <- layout.l2[,2] + (l)*layer.shift.y
        
        z1 <- (l-1)*layer.space
        z2 <- (l)*layer.space
        
        # Identificar nodos activos en ambas
        act1 <- which(igraph::degree(g.list[[l]]) > 0)
        act2 <- which(igraph::degree(g.list[[l+1]]) > 0)
        common <- intersect(act1, act2)
        
        if(length(common)>0){
             x.from <- layout.l1[common,1]
             y.from <- layout.l1[common,2]
             x.to   <- layout.l2[common,1]
             y.to   <- layout.l2[common,2]
             
             rgl::segments3d(as.vector(rbind(x.from, x.to)), as.vector(rbind(y.from, y.to)), as.vector(rbind(rep(z1, length(common)), rep(z2, length(common)))),
                             col="gray80", alpha=0.3)
        }
    }
}


# --- 4. VISUALIZACIÓN 3D CON MUXVIZ (REAL) ---
if (requireNamespace("rgl", quietly = TRUE)) {
  library(rgl)
  # Necesario para que RMarkdown renderice WebGL
  setupKnitr(autoprint = TRUE)
  
  # Configuración de colores y layout
  colores <- c(bicimad="#1a9641", bus="#d7191c", metro="#2c7bb6", parking="#fdae61")
  
  # Asegurar que el layout es una matriz numérica Nx2
  # IMPORTANTE: Normalizamos las coordenadas para que la escala Z (spacing=2) sea visible
  # Las coordenadas UTM son enormes (400000, 4400000), haciendo que Z=2 sea invisible.
  layout_raw <- as.matrix(centroides_mat)
  layout_geo <- scale(layout_raw) # Centra y escala (aprox -2 a 2)
  
  # MuxViz nativo falla en R 4.3+ por verificación de tipos en vectores (length > 1)
  # Usamos nuestra versión corregida basada en muxViz:
  
  plot_multiplex3D_fix(
    g_list, 
    layer.layout = layout_geo, 
    layer.colors = colores,
    layer.shift.x = 0,   # Sin shift horizontal extra
    layer.shift.y = 0,
    layer.space = 4,
  )
  
  # Ajustar vista
  view3d(theta = 45, phi = 25, zoom = 0.6)
  
  # Renderizar widget interactivo
  rglwidget()
  
} else {
  message("Instala rgl para ver el gráfico 3D")
}

Visualización de la red por capas y agregada (2D):

# --- CONFIGURACIÓN ESTÉTICA ---
colores <- c(bicimad="#1a9641", bus="#d7191c", metro="#2c7bb6", parking="#fdae61")
titulos <- names(g_list)
# Layout geográfico fijo (usamos las coordenadas que ya calculaste)
layout_geo <- as.matrix(centroides_mat)


# ==============================================================================
# VISTA 1: PANEL 2D (MULTIPLEX SLICES)
# Muestra cada capa por separado manteniendo la posición geográfica
# ==============================================================================
par(mfrow = c(2, 2), mar = c(1, 1, 3, 1)) # Rejilla 2x2

for (m in titulos) {
  g <- g_list[[m]]
  deg <- degree(g)
  
  # Lógica visual: 
  # Nodos activos = color intenso y tamaño según grado
  # Nodos inactivos = gris muy claro y pequeños (para mantener referencia espacial)
  v_col <- ifelse(deg > 0, colores[m], adjustcolor("gray90", alpha.f = 0.5))
  v_size <- ifelse(deg > 0, 3 + (deg/max(deg))*5, 1)
  v_frame <- ifelse(deg > 0, "black", NA) # Sin borde si está inactivo
  
  # Aristas: Más finas y transparentes
  e_col <- adjustcolor(colores[m], alpha.f = 0.6)
  
  plot(g,
       layout = layout_geo,
       vertex.label = NA,
       vertex.size = v_size,
       vertex.color = v_col,
       vertex.frame.color = v_frame,
       edge.color = e_col,
       edge.width = 0.8,
       edge.arrow.size = 0,
       main = toupper(m))
  
  box(col = "gray80") # Marco suave alrededor de cada plot
}

# ==============================================================================
# VISTA 2: RED AGREGADA (OVERLAY)
# Superpone todo para ver "corredores" multimodales
# ==============================================================================
par(mfrow = c(1, 1), mar = c(0, 0, 2, 0))

# Creamos un plot vacío primero para establecer los límites
plot(layout_geo, type = "n", axes = FALSE, xlab = "", ylab = "", 
     main = "RED AGREGADA (OVERLAP)")

# Dibujamos las aristas de cada capa acumulativamente
# Truco: Dibujamos primero las capas más densas (Bus/Metro) al fondo, Bicimad arriba
orden_capas <- c("bus", "metro", "parking", "bicimad") 

for (m in orden_capas) {
  if (m %in% names(g_list)) {
    g <- g_list[[m]]
    
    # Extraemos aristas solo si existen
    if (gsize(g) > 0) {
      el <- as_edgelist(g, names = FALSE)
      
      # Coordenadas de inicio y fin de las aristas
      x0 <- layout_geo[el[,1], 1]
      y0 <- layout_geo[el[,1], 2]
      x1 <- layout_geo[el[,2], 1]
      y1 <- layout_geo[el[,2], 2]
      
      # Dibujar líneas (curvadas ligeramente para ver solapamientos sería ideal, 
      # pero rectas es más limpio en mapas densos)
      segments(x0, y0, x1, y1, 
               col = adjustcolor(colores[m], alpha.f = 0.4), 
               lwd = 1)
    }
  }
}

# Dibujamos nodos ENCIMA de las líneas
# Calculamos grado total sumando todas las capas
deg_matrix <- sapply(g_list, degree)
deg_total <- rowSums(deg_matrix)
active_nodes <- deg_total > 0

points(layout_geo[active_nodes, ], 
       pch = 21, 
       bg = "white", 
       col = "gray30", 
       cex = 0.8 + (deg_total[active_nodes]/max(deg_total))*1.5)

legend("topright", legend = names(colores), 
       col = colores, lty = 1, lwd = 3, bty = "n", cex = 0.8)

Finalmente, para poder tener una perspectiva distinta utilizamos el resto de la visualización en un entorno tridimensional interactivo. A diferencia de las proyecciones planas, este gráfico ‘desapila’ la red de transporte ubicando cada modo (Metro, Bus, BiciMAD, Parking) en un nivel diferente del eje Z. Esto nos permite observar la estructura multiplex del sistema: las líneas verticales negras conectan una misma zona geográfica a través de sus diferentes capas, facilitando la identificación visual de los nodos con mayor integración multimodal.

# ===========================================
# Visualización 3D multicapa (inter-capa con scatter3d lines)
# ===========================================

if (requireNamespace("plotly", quietly = TRUE)) {

  library(plotly)
  library(scales)

  # ---- 1. Capas y eje Z
  layer_order <- modos
  layer_index <- setNames(seq_along(layer_order) - 1, layer_order)  # 0,1,2,3

  # ---- 2. Nodos multilayer: solo zonas donde ese modo existe
  nodes_list <- list()
  for (m in layer_order) {
    # Filtramos zonas con ese modo activo
    sub_zonas <- zonas_attrs[zonas_attrs[[m]] == 1, ]
    
    # Cruzamos con centroides
    df_m <- merge(sub_zonas[, "cluster_id", drop = FALSE], zonas_centroides, by = "cluster_id")
    
    # Construimos el dataframe de la capa
    capa_df <- data.frame(
      node       = paste(df_m$cluster_id, m, sep = "_"),
      cluster_id = df_m$cluster_id,
      layer      = m,
      x_raw      = df_m$x,
      y_raw      = df_m$y,
      stringsAsFactors = FALSE
    )
    nodes_list[[m]] <- capa_df
  }
  nodes_multi_3d <- do.call(rbind, nodes_list)

  # Normalizamos X, Y y asignamos Z
  nodes_multi_3d$x <- rescale(nodes_multi_3d$x_raw, to = c(0, 1))
  nodes_multi_3d$y <- rescale(nodes_multi_3d$y_raw, to = c(0, 1))
  nodes_multi_3d$z <- layer_index[nodes_multi_3d$layer]

  # ---- 3. Multimodalidad por zona
  # Calculamos cuántos modos hay activos sumando las columnas de modos
  modos_activos_vec <- rowSums(zonas_attrs[, modos])
  node_multimodalidad <- data.frame(
    cluster_id = zonas_attrs$cluster_id,
    modos_activos = modos_activos_vec
  )

  # Clusters con al menos 2 modos
  clusters_interes <- node_multimodalidad$cluster_id[node_multimodalidad$modos_activos >= 2]

  # ---- 4. Aristas inter-capa: todos los pares de modos presentes
  edges_inter_list <- list()
  for (cid in clusters_interes) {
    fila <- zonas_attrs[zonas_attrs$cluster_id == cid, modos]
    modos_presentes <- modos[as.logical(unlist(fila))]

    if (length(modos_presentes) >= 2) {
      comb <- t(combn(modos_presentes, 2))
      edges_inter_list[[as.character(cid)]] <- data.frame(
        from       = paste(cid, comb[, 1], sep = "_"),
        to         = paste(cid, comb[, 2], sep = "_"),
        cluster_id = cid,
        stringsAsFactors = FALSE
      )
    }
  }
  edges_inter_3d <- do.call(rbind, edges_inter_list)

  # ---- 5. Posiciones de nodos para las aristas
  pos_nodes <- nodes_multi_3d[, c("node", "x", "y", "z", "layer")]

  edges_plot <- merge(edges_inter_3d, pos_nodes, by.x = "from", by.y = "node")
  edges_plot <- merge(edges_plot, pos_nodes, by.x = "to", by.y = "node", suffixes = c("", "_to"))
  
  # Renombrar manualmente para claridad
  names(edges_plot)[names(edges_plot) == "x_to"] <- "xend"
  names(edges_plot)[names(edges_plot) == "y_to"] <- "yend"
  names(edges_plot)[names(edges_plot) == "z_to"] <- "zend"

  # ---- 6. Convertimos aristas a formato "long" para scatter3d lines
  if (nrow(edges_plot) > 0) {
    # Creamos vectores con NAs intermedios para separar líneas en plotly
    x_long <- as.vector(t(cbind(edges_plot$x, edges_plot$xend, NA)))
    y_long <- as.vector(t(cbind(edges_plot$y, edges_plot$yend, NA)))
    z_long <- as.vector(t(cbind(edges_plot$z, edges_plot$zend, NA)))
    edges_long <- data.frame(x = x_long, y = y_long, z = z_long)
  } else {
    edges_long <- data.frame(x = numeric(), y = numeric(), z = numeric())
  }


  # ---- 7. Añadimos multimodalidad y tamaño de nodo
  nodes_multi_3d <- merge(nodes_multi_3d, node_multimodalidad, by = "cluster_id", all.x = TRUE)
  nodes_multi_3d$modos_activos[is.na(nodes_multi_3d$modos_activos)] <- 1
  nodes_multi_3d$size <- rescale(nodes_multi_3d$modos_activos, to = c(4, 11))

  # ---- 8. Colores por capa
  layer_colors <- c(
    bicimad = "#1a9641",
    bus     = "#d7191c",
    metro   = "#2c7bb6",
    parking = "#fdae61"
  )

  # ---- 9. Gráfico 3D
  p3d <- plot_ly()

  # Nodos
  p3d <- add_trace(
    p3d,
    data = nodes_multi_3d,
    type = "scatter3d",
    mode = "markers",
    x = ~x, y = ~y, z = ~z,
    color = ~layer,
    colors = layer_colors,
    marker = list(size = ~size),
    hoverinfo = "text",
    text = ~paste(
      "Zona:", cluster_id,
      "<br>Modo:", layer,
      "<br>Modos activos (total cluster):", modos_activos
    )
  )

  # Aristas inter-capa como líneas negras
  if (nrow(edges_long) > 0) {
    p3d <- add_trace(
      p3d,
      data = edges_long,
      type = "scatter3d",
      mode = "lines",
      opacity = 0.1,
      x = ~x, y = ~y, z = ~z,
      line = list(color = "darkgray", width = 4),
      showlegend = FALSE,
      hoverinfo = "none"
    )
  }

  p3d <- layout(
    p3d,
    title = "Red multicapa de transporte de Madrid (3D, enlaces entre capas)",
    scene = list(
      xaxis = list(title = "X (normalizada)"),
      yaxis = list(title = "Y (normalizada)"),
      zaxis = list(
        title = "Capa (modo de transporte)",
        tickvals = layer_index,
        ticktext = layer_order
      ),
      aspectmode = "cube"
    )
  )

  p3d

} else {
  message("Instala 'plotly'")
}

Interpretación de la visualización 3D multicapa

La visualización tridimensional obtenida permite observar de forma explícita la estructura multicapa del sistema de transporte urbano de Madrid, donde cada plano horizontal representa un modo de transporte (bicimad, bus, metro y parking) y las conexiones verticales indican zonas de intercambio multimodal.

Separación clara de capas y estructura espacial

En primer lugar, se aprecia una separación nítida entre las cuatro capas, lo que confirma que cada modo de transporte presenta una distribución espacial propia:

La capa bicimad aparece más concentrada espacialmente, reflejando su carácter urbano y centralizado.
La capa bus muestra una cobertura más extensa, coherente con su función de red capilar que conecta barrios periféricos.
La capa metro se sitúa entre ambas, con una estructura más compacta que refleja la jerarquía de estaciones y líneas.
La capa parking aparece claramente diferenciada y menos densa, lo que concuerda con su función de infraestructura de apoyo más localizada.

Esta separación valida la decisión metodológica de modelar el sistema como una red multicapa, ya que una red agregada ocultaría estas diferencias estructurales.

Uniones verticales como indicador de multimodalidad

Las líneas verticales que conectan nodos entre capas representan zonas donde coexisten varios modos de transporte. Estas conexiones no aparecen de forma homogénea, sino que se concentran en un subconjunto reducido de zonas, lo que pone de manifiesto una fuerte heterogeneidad en la multimodalidad del sistema.

En términos de teoría de redes:

Estas zonas actúan como nodos altamente versátiles, ya que participan simultáneamente en varias capas.
Son puntos críticos para la integración funcional del sistema, permitiendo el transbordo eficiente entre modos.
Su posición estructural las convierte en puntos de control del flujo global de movilidad.

Visualmente, estas zonas destacan como “columnas” que atraviesan varias capas, reforzando la idea de que la conectividad global del sistema depende de un número limitado de nodos clave.

Versatilidad y roles multinodales

La figura ilustra de manera intuitiva el concepto de versatilidad de nodo:

La mayoría de los nodos aparecen en una sola capa, actuando como nodos especializados (por ejemplo, paradas de bus o estaciones de bici aisladas).
En contraste, un número reducido de nodos conecta tres o incluso cuatro capas, desempeñando un rol multinodal esencial.

Este patrón es consistente con la literatura sobre redes multicapas, donde se observa que la importancia de un nodo no puede evaluarse correctamente desde una sola capa. La visualización 3D confirma que los nodos versátiles no solo existen, sino que son estructuralmente visibles cuando se representan explícitamente las capas y sus interdependencias.

Implicaciones para vulnerabilidad y resiliencia

Desde el punto de vista de la vulnerabilidad del sistema, la visualización sugiere una conclusión relevante:

La conectividad multimodal del transporte madrileño depende de un conjunto relativamente pequeño de nodos altamente interconectados entre capas.
Una perturbación localizada en estos nodos (por ejemplo, cierre de un intercambiador o fallo simultáneo de varios modos en una zona) podría propagarse rápidamente al resto del sistema, afectando a múltiples capas.

Este comportamiento es coherente con los modelos de cascadas de fallo en redes interdependientes descritos en la literatura (Buldyrev et al., 2010; Duan et al., 2019), donde la interdependencia aumenta tanto la eficiencia como la fragilidad del sistema.

Valor añadido de la representación multicapa

Finalmente, esta visualización demuestra el valor añadido del enfoque multicapa frente a representaciones tradicionales:

Permite identificar de forma directa los intercambiadores críticos.
Hace visible la estructura de interdependencia entre modos.
Facilita la interpretación de conceptos abstractos como versatilidad, roles multinodales y vulnerabilidad sistémica.

En conjunto, el gráfico 3D no solo actúa como herramienta exploratoria, sino como una evidencia visual clara de que el sistema de transporte de Madrid funciona como una red interconectada de capas, cuya eficiencia y resiliencia dependen de la correcta gestión de sus nodos multimodales clave.

5 Centralidad y versatilidad de las zonas

5.1 Medidas de centralidad por capa

Calculamos grado y fuerza (suma de pesos) para cada zona en cada modo.

# --- 1. PREPARACIÓN DE DATOS PARA MUXVIZ ---


# Mapeamos Modos a IDs numéricos (1..L)
map_modos <- data.frame(modo = modos, layer_id = 1:length(modos))

# Mapeamos Zonas (Cluster IDs) a IDs de nodo consecutivos (1..N)
# Esto es crucial porque muxViz espera índices matriciales 1..N
mis_nodos <- sort(unique(zonas_attrs$cluster_id))
map_nodos <- data.frame(
  cluster_id = mis_nodos,
  node_id    = 1:length(mis_nodos)
)

# Creamos la "Extended Edge List" para aristas INTRA-CAPA
edges_muxviz_intra <- edges_intra %>%
  left_join(map_modos, by = "modo") %>%
  rename(layer_from = layer_id) %>%
  mutate(layer_to = layer_from) %>%
  left_join(map_nodos, by = c("from" = "cluster_id")) %>%
  rename(node_from = node_id) %>%
  left_join(map_nodos, by = c("to" = "cluster_id")) %>%
  rename(node_to = node_id) %>%
  select(node_from, layer_from, node_to, layer_to, w)

# Creamos aristas INTER-CAPA (Coupling) de forma vectorizada (mucho más rápido)
nodos_capas_activas <- zonas_attrs %>%
  select(cluster_id, all_of(modos)) %>%
  pivot_longer(cols = all_of(modos), names_to = "modo", values_to = "active") %>%
  filter(active == 1) %>%
  left_join(map_modos, by = "modo") %>%
  left_join(map_nodos, by = "cluster_id") %>%
  select(node_id, layer_id)

edges_muxviz_inter <- nodos_capas_activas %>%
  rename(layer_from = layer_id) %>%
  inner_join(nodos_capas_activas %>% rename(layer_to = layer_id), 
             by = "node_id", relationship = "many-to-many") %>%
  filter(layer_from < layer_to) %>% # Solo pares únicos
  mutate(
    node_from = node_id,
    node_to = node_id,
    w = 1 # Peso estándar para cambio de modo
  ) %>%
  select(node_from, layer_from, node_to, layer_to, w)

# Combinamos y forzamos data.frame
edges_full <- as.data.frame(rbind(edges_muxviz_intra, edges_muxviz_inter))

# --- 2. CONSTRUCCIÓN DE LA SUPRA-MATRIZ (MANUAL ROBUSTA) ---
Layers <- length(modos)
Nodes  <- nrow(map_nodos)
TotalNodes <- Layers * Nodes

get_supra_idx <- function(n, l, NumNodes) {
  return( (l - 1) * NumNodes + n )
}

edges_supra <- edges_full %>%
  mutate(
    idx_from = get_supra_idx(node_from, layer_from, Nodes),
    idx_to   = get_supra_idx(node_to, layer_to, Nodes)
  )

# Construimos matriz esparsa
library(Matrix)
SupraAdj <- sparseMatrix(
  i = edges_supra$idx_from,
  j = edges_supra$idx_to,
  x = edges_supra$w,
  dims = c(TotalNodes, TotalNodes),
  symmetric = FALSE
)

# Aseguramos simetría
SupraAdj <- SupraAdj + t(SupraAdj)

# Conversión final a formato igraph adjacency para asegurar compatibilidad total
g_supra_temp <- graph_from_data_frame(edges_supra[, c("idx_from", "idx_to", "w")], directed = FALSE, vertices = 1:TotalNodes)
SupraAdj <- as_adjacency_matrix(g_supra_temp, attr = "w", sparse = TRUE)

# Limpieza de memoria agresiva (objetos grandes ya no necesarios)
rm(edges_supra, edges_muxviz_intra, edges_muxviz_inter, g_supra_temp)
gc()

##           used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells 3057524 163.3   10232512 546.5         NA  7962295 425.3
## Vcells 6476556  49.5   15157805 115.7      16384 15157786 115.7

# --- 3. CÁLCULO DE MÉTRICAS MULTICAPA (MUXVIZ) ---
# Grado Multicapa
# Grado (Robusto)
tryCatch({
  multi_degree <- GetMultiDegree(SupraAdj, Layers = Layers, Nodes = Nodes, isDirected = FALSE)
}, error = function(e) multi_degree <<- rowSums(SupraAdj))

# PageRank (Robusto con igraph sobre supra-matriz)
# Usamos igraph::page_rank que es C++ optimizado (instantáneo y robusto)
g_supra_final <- graph_from_adjacency_matrix(SupraAdj, mode = "undirected", weighted = TRUE)
pr_result <- page_rank(g_supra_final, damping = 0.85)

# Agregamos los valores de las réplicas para obtener el valor del nodo físico
# El vector 'pr_result$vector' tiene longitud Layers * Nodes
pr_vals <- pr_result$vector
node_indices <- rep(1:Nodes, Layers) # Vector de índices de nodo repetidos
multi_pagerank <- tapply(pr_vals, node_indices, sum) # Suma por nodo físico

# Consolidamos resultados en un DF
metricas_muxviz <- data.frame(
  node_id = 1:Nodes,
  multi_degree_muxviz = as.numeric(multi_degree),
  multi_pagerank_muxviz = as.numeric(multi_pagerank)
) %>%
  left_join(map_nodos, by = "node_id")

# --- 4. INTEGRACIÓN CON DATOS ANTERIORES ---
centralidad_lista <- list()
for (modo in names(g_list)) {
  g <- g_list[[modo]]
  centralidad_lista[[modo]] <- data.frame(
    cluster_id = as.integer(V(g)$name),
    modo       = modo,
    degree     = degree(g), # Grado local (intra-capa)
    strength   = strength(g, weights = E(g)$w),
    stringsAsFactors = FALSE
  )
}
centralidad_por_capa <- do.call(rbind, centralidad_lista)

kable(head(centralidad_por_capa), caption = "Centralidad intra-capa (calculada con igraph)")

Centralidad intra-capa (calculada con igraph)
	cluster_id	modo	degree	strength
bicimad.1	1	bicimad	5	0.0071526
bicimad.2	2	bicimad	6	0.0131310
bicimad.3	3	bicimad	4	0.0072934
bicimad.4	4	bicimad	6	0.0105906
bicimad.5	5	bicimad	7	0.0130558
bicimad.6	6	bicimad	6	0.0151738

Distribución del grado por modo (intra-capa):

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
for (m in modos) {
  datos_sub <- centralidad_por_capa[centralidad_por_capa$modo == m, ]
  hist(datos_sub$degree, 
       breaks = 30, 
       col = adjustcolor(colores[m], alpha.f = 0.6),
       border = "white",
       main = paste("MODO:", toupper(m)),
       xlab = "Grado Intra-capa", ylab = "Frecuencia")
}

5.2 Versatilidad y Ranking Multicapa (con muxViz)

Utilizando la librería muxViz, calculamos métricas globales que integran la estructura topológica de todas las capas simultáneamente.

Versatilidad de Grado (Entropía): Mide el equilibrio de conexiones entre capas.
Multi-Layer PageRank: Mide la importancia global del nodo considerando flujos que viajan a través de múltiples modos.

# A. Versatilidad Manual (Entropía)
mat_grado_df <- as.data.frame.matrix(xtabs(degree ~ cluster_id + modo, data = centralidad_por_capa))
for (m in modos) if (!m %in% names(mat_grado_df)) mat_grado_df[[m]] <- 0
grado_mat <- as.matrix(mat_grado_df[, modos])

calc_versatilidad <- function(vec) {
  if (sum(vec) == 0) return(NA_real_)
  p <- vec / sum(vec)
  p <- p[p > 0]
  H <- -sum(p * log(p))
  H / log(length(vec))
}
versatilidad_entropia <- apply(grado_mat, 1, calc_versatilidad)

# B. Unificación en Dataframe Final
# Unimos la entropía manual con el PageRank calculado por muxViz
versatilidad_df <- data.frame(
  cluster_id = as.integer(rownames(grado_mat)),
  versatilidad_grado = as.numeric(versatilidad_entropia)
) %>%
  left_join(metricas_muxviz, by = "cluster_id") %>% # Añadimos PageRank de muxViz
  left_join(zonas_attrs, by = "cluster_id") %>%
  left_join(viajes_distrito, by = c("distrito_asignado" = "distrito"))

summary(versatilidad_df$multi_pagerank_muxviz)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0002212 0.0002974 0.0003952 0.0006418 0.0006399 0.0027629

5.2.1 Nodos más versátiles y más especializados

# --- 1. TOP 15 ZONAS MÁS VERSÁTILES ---
orden_desc <- order(versatilidad_df$versatilidad_grado, decreasing = TRUE)
top_versatiles_full <- head(versatilidad_df[orden_desc, ], n = 15)

# Verificamos qué columnas están disponibles para evitar el error "undefined columns"
cols_deseadas <- c("cluster_id", "versatilidad_grado", "modos_disponibles", "distrito_asignado")
cols_presentes <- cols_deseadas[cols_deseadas %in% colnames(top_versatiles_full)]

top_versatiles <- top_versatiles_full[, cols_presentes]

# --- 2. TOP 15 ZONAS MÁS ESPECIALIZADAS ---
orden_asc <- order(versatilidad_df$versatilidad_grado, na.last = TRUE)
top_especializados_full <- head(versatilidad_df[orden_asc, ], n = 15)

top_especializados <- top_especializados_full[, cols_presentes]

# --- 3. TABLAS RESULTADO ---
kable(top_versatiles, caption = "Top 15 zonas más versátiles")

Top 15 zonas más versátiles
	cluster_id	versatilidad_grado	distrito_asignado
164	164	0.9983119	2807907
446	446	0.9976183	2807904
101	101	0.9968869	2807907
41	41	0.9955380	2807907
56	56	0.9949672	2807904
283	283	0.9949672	2807903
239	239	0.9946468	2807901
814	814	0.9934804	2807904
71	71	0.9927376	2807904
339	339	0.9927376	2807911
534	534	0.9927376	2807908
7	7	0.9926388	2807901
251	251	0.9892423	2807916
562	562	0.9892423	2807920
573	573	0.9863237	2807904

kable(top_especializados, caption = "Top 15 zonas más especializadas")

Top 15 zonas más especializadas
	cluster_id	distrito_asignado
9	9	2807903
10	10	2807909
15	15	2807921
17	17	2807921
19	19	2807921
20	20	2807921
21	21	2807921
22	22	2807908
26	26	2807911
31	31	2807921
32	32	2807921
33	33	2807921
34	34	2807911
36	36	2807911
38	38	2807920

Desde el punto de vista espacial, las zonas con alta versatilidad de grado se concentran en torno a los grandes intercambiadores de transporte y a las áreas más centrales de la ciudad. Se trata de nodos que combinan conectividad en varias capas de la red (metro, bus, cercanías, bicicleta pública), de modo que actúan como puntos de transferencia entre modos diferentes. Estas zonas tienen grados relativamente altos en varias capas simultáneamente y, por tanto, desempeñan un papel de “puente multimodal” entre distintas partes del sistema de transporte.

En el extremo opuesto, las zonas con baja versatilidad tienden a estar especializadas en un único modo de transporte, típicamente la red de autobuses o una línea concreta de metro. Su grado puede ser elevado dentro de una sola capa, pero su conectividad se limita a ese modo, de modo que su capacidad para redistribuir flujos entre capas es reducida. Esta especialización puede ser eficiente en términos operativos, pero hace que la zona sea menos flexible ante perturbaciones que afecten a su modo dominante, ya que existen menos alternativas multimodales inmediatas.

5.3 Relación entre versatilidad e intensidad de viajes

# Quitamos NAs de las columnas que vamos a usar
datos_plot <- na.omit(versatilidad_df[, c("viajes_medios_por_persona", "versatilidad_grado")])

# Ordenamos por el eje X (necesario para que la línea no se cruce)
datos_plot <- datos_plot[order(datos_plot$viajes_medios_por_persona), ]

# Gráfico base
plot(datos_plot$viajes_medios_por_persona, datos_plot$versatilidad_grado,
     pch = 19, col = adjustcolor("black", 0.3),
     xlab = "Viajes medios por persona", ylab = "Versatilidad")

# Línea de tendencia (loess)
modelo <- loess(versatilidad_grado ~ viajes_medios_por_persona, data = datos_plot)
lines(datos_plot$viajes_medios_por_persona, predict(modelo), col = "blue", lwd = 2)

La Figura anterior muestra la relación entre la y el a nivel de distrito. Visualmente, la nube de puntos no presenta una tendencia clara, lo que sugiere una relación débil entre ambas variables. Este resultado se confirma mediante el análisis de correlación de Pearson.

datos_cor <- versatilidad_df %>% 
  filter(!is.na(versatilidad_grado),
         !is.na(viajes_medios_por_persona))

n_obs <- nrow(datos_cor)
n_obs

## [1] 1558

if (n_obs >= 3) {
  cor_test <- cor.test(
    datos_cor$versatilidad_grado,
    datos_cor$viajes_medios_por_persona
  )
  cor_test
} else {
  cat("No hay suficientes observaciones completas para estimar de forma fiable la correlación entre versatilidad y viajes.
")
}

## 
##  Pearson's product-moment correlation
## 
## data:  datos_cor$versatilidad_grado and datos_cor$viajes_medios_por_persona
## t = 0.14253, df = 1556, p-value = 0.8867
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.04605711  0.05326588
## sample estimates:
##         cor 
## 0.003613297

Este cálculo arroja un coeficiente prácticamente nulo (\(r = 0.0036\)) y no significativo (\(p = 0.8867\), \(n = 1558\)). El intervalo de confianza al 95,% incluye el valor cero, lo que indica la ausencia de una asociación lineal estadísticamente significativa entre ambas magnitudes.

La línea de suavizado refuerza esta interpretación, mostrando un comportamiento prácticamente plano a lo largo del rango de valores observados. Aunque se aprecia una ligera curvatura en los valores más altos de viajes medios, esta variación es reducida y carece de relevancia estadística. En consecuencia, no se observa evidencia de que los distritos con mayor intensidad de movilidad presenten zonas con una conectividad más equilibrada entre modos de transporte.

Este resultado es coherente con la definición de la versatilidad de grado como una medida basada en la distribución del grado entre capas modales, y no en el volumen absoluto de actividad. Una zona puede registrar un número elevado de viajes concentrados en un único modo y, por tanto, mostrar una baja versatilidad, mientras que otras zonas con menor intensidad de uso pueden desempeñar un papel multimodal más destacado. De este modo, la versatilidad de grado captura una dimensión complementaria de la red, asociada a su diversidad modal, que no depende directamente de la cantidad de viajes generados.

6 Vulnerabilidad de la red de transporte

En lugar de implementar un modelo muy complejo de cascadas, realizamos un análisis sencillo pero informativo de robustez estructural:

Tomamos la red agregada (unión de todas las capas).
Eliminamos nodos de forma dirigida (los más centrales) en una capa y observamos cómo disminuye el tamaño de la componente gigante agregada.

6.1 Red agregada

# Matrices de adyacencia por capa
adj_list <- lapply(g_list, function(g) {
  as.matrix(as_adj(g, attr = "w", sparse = FALSE))
})

# Aseguramos mismo orden de nodos
orden_nodos <- order(as.integer(V(g_list[[1]])$name))
adj_list <- lapply(adj_list, function(m) m[orden_nodos, orden_nodos])

adj_agregada <- Reduce("+", adj_list)

g_agregado <- graph_from_adjacency_matrix(adj_agregada, mode = "undirected", weighted = TRUE)

comp <- components(g_agregado)
tamano_gcc_inicial <- max(comp$csize)
tamano_gcc_inicial

## [1] 1528

6.2 Robustez frente a fallos en una capa

Analizamos qué ocurre si vamos eliminando nodos con mayor grado en la capa metro y miramos el tamaño relativo de la componente gigante en la red agregada.

# --- 1. PREPARACIÓN DEL ORDEN DE FALLO ---
g_metro <- g_list[["metro"]]
deg_metro <- degree(g_metro)

# Nombres de los nodos ordenados de mayor a menor grado
orden_fallo <- names(sort(deg_metro, decreasing = TRUE))

# --- 2. SIMULACIÓN DE ROBUSTEZ ---
fracciones <- seq(0, 0.5, by = 0.05)
res_lista <- list()

for (i in seq_along(fracciones)) {
  f <- fracciones[i]
  
  # Calculamos cuántos nodos eliminar
  n_remove <- ceiling(f * length(orden_fallo))
  nodos_fuera <- if (n_remove > 0) orden_fallo[1:n_remove] else character(0)
  
  # Eliminamos los nodos de la red agregada y calculamos el GCC
  g_agregado_reducido <- delete_vertices(g_agregado, nodos_fuera)
  comp_red <- components(g_agregado_reducido)
  
  gcc_red <- if (length(comp_red$csize) == 0) 0 else max(comp_red$csize)
  
  # Guardamos resultados en un dataframe temporal
  res_lista[[i]] <- data.frame(
    frac_eliminada = f,
    tamano_gcc_rel = gcc_red / tamano_gcc_inicial
  )
}

# Unimos los resultados
resultado_robustez <- do.call(rbind, res_lista)

# --- 3. VISUALIZACIÓN ---
# Configuramos el gráfico base
plot(resultado_robustez$frac_eliminada, resultado_robustez$tamano_gcc_rel,
     type = "b", # Tipo "both" (puntos y líneas)
     pch = 19, 
     col = "#2c7bb6", # Azul Metro
     lwd = 2,
     main = "Robustez ante fallos en la capa METRO",
     xlab = "Fracción de nodos eliminados (%)",
     ylab = "Tamaño relativo del GCC (%)",
     xaxt = "n", yaxt = "n",
     frame.plot = FALSE)

# Personalizamos los ejes para mostrar porcentajes
eje_x <- seq(0, 0.5, 0.1)
axis(1, at = eje_x, labels = paste0(eje_x * 100, "%"))

eje_y <- seq(0, 1, 0.2)
axis(2, at = eje_y, labels = paste0(eje_y * 100, "%"), las = 1)

grid(col = "gray90", lty = "solid")

La Figura muestra la evolución del tamaño relativo de la componente gigante de la red agregada a medida que se eliminan progresivamente nodos de la capa de metro, siguiendo un criterio dirigido basado en su centralidad. Se observa una disminución aproximadamente monótona y casi lineal de la conectividad global, lo que indica que los nodos más centrales de la red de metro desempeñan un papel estructural clave en la cohesión del sistema multicapa.

En particular, la eliminación del 20–30 % de los nodos más relevantes de la capa de metro reduce el tamaño de la componente gigante a alrededor del 70–75 % de su valor inicial, mientras que al eliminar la mitad de los nodos centrales la red pierde más del 50 % de su conectividad global. Este resultado pone de manifiesto que la capa de metro actúa como columna vertebral del sistema de transporte urbano: disrupciones severas concentradas en sus nodos más importantes pueden propagarse al resto de capas y provocar una fragmentación significativa de la red agregada.

6.2.1 Comparación entre fallos dirigidos y aleatorios

set.seed(123)

simular_ataque <- function(g, estrategia = c("dirigido", "aleatorio"),
                           fracciones = seq(0, 0.5, by = 0.05)) {
  estrategia <- match.arg(estrategia)
  n <- vcount(g)
  resultados <- tibble()
  
  if (estrategia == "dirigido") {
    orden <- names(sort(degree(g), decreasing = TRUE))
  } else {
    orden <- sample(V(g)$name)
  }
  
  for (f in fracciones) {
    k <- ceiling(f * n)
    g_tmp <- delete_vertices(g, orden[seq_len(k)])
    gcc <- if (vcount(g_tmp) == 0) 0 else max(components(g_tmp)$csize)
    
    resultados <- bind_rows(
      resultados,
      tibble(
        frac_eliminada = f,
        tamano_gcc_rel = gcc / tamano_gcc_inicial,
        estrategia = estrategia
      )
    )
  }
  resultados
}

res_dirigido  <- simular_ataque(g_agregado, "dirigido")
res_aleatorio <- simular_ataque(g_agregado, "aleatorio")

res_vuln <- bind_rows(res_dirigido, res_aleatorio)

ggplot(res_vuln,
       aes(x = frac_eliminada, y = tamano_gcc_rel, color = estrategia)) +
  geom_line(size = 1) +
  geom_point() +
  theme_minimal() +
  labs(
    x = "Fracción de nodos eliminados",
    y = "Tamaño relativo de la componente gigante",
    color = "Estrategia",
    title = "Robustez de la red: ataque dirigido vs fallo aleatorio"
  )

El análisis conjunto de la versatilidad nodal y la robustez estructural de la red agregada de Madrid arroja tres conclusiones fundamentales sobre la arquitectura del sistema de transporte:

Independencia entre Infraestructura y Demanda Individual. Tal como se observa en el gráfico de Versatilidad vs. Viajes, no existe una correlación clara entre la intensidad de uso (viajes medios por persona en el distrito) y la complejidad multimodal de las paradas.
- Estratificación de la oferta: La disposición de los puntos en bandas horizontales discretas revela que la red no es continua, sino jerárquica: la gran mayoría de nodos son unimodales (base del gráfico), mientras que la multimodalidad se da “a saltos” (agregando capas de Metro o BiciMAD).
- Homogeneidad de la demanda: El rango estrecho del eje X (1.95 - 2.15 viajes) confirma que la movilidad individual es uniforme en toda la ciudad. Por tanto, la ubicación de los grandes intercambiadores (nodos con alta versatilidad) responde a criterios estratégicos de planificación de red y transbordo, y no a una mayor intensidad de viajes per cápita en esos distritos específicos.
La red es “Robusta pero Frágil”. La comparación entre estrategias de fallo (gráfico de Ataque dirigido vs. Fallo aleatorio) confirma que la red de transporte de Madrid comparte la topología típica de las redes complejas libres de escala (scale-free):
- Alta tolerancia a fallos aleatorios (Línea Roja): El sistema es extremadamente resiliente ante incidencias fortuitas. Incluso ante la inoperatividad del 40% de los nodos de forma aleatoria (simulando averías dispersas o cortes menores), la red mantiene una componente gigante conectada significativa, garantizando la movilidad en gran parte del territorio.
- Extrema fragilidad ante ataques dirigidos (Línea Azul/Turquesa): Por el contrario, la red es muy vulnerable si los fallos afectan a sus nodos centrales. La eliminación de apenas el 15-20% de los nodos más conectados provoca el colapso funcional del sistema (fragmentación de la componente gigante por debajo del 20%), aislando grandes zonas de la ciudad.
El Metro como columna vertebral crítica. El análisis específico de la capa de Metro (gráfico de caída monótona) demuestra su rol de esqueleto estructural. Al eliminar los nodos principales de esta capa, la conectividad global de la red agregada decae casi linealmente. Esto implica que la cohesión entre los distintos modos (bus, bici, parking) depende críticamente de la operatividad de las estaciones de metro centrales, que actúan como los “pegamentos” que unen las distintas capas de la red multicapa.

7 Dashboard Interactivo de la Red de Transporte

En esta sección se presenta un dashboard interactivo que resume las métricas principales del análisis de redes multicapas de Madrid. Este dashboard permite una exploración visual e interactiva de los resultados obtenidos.

7.1 Métricas Generales del Sistema

7.2 Distribución de Nodos por Capa

Análisis de la Distribución: La preponderancia de la capa de Autobús (representando la mayor fracción de nodos activos) confirma su rol de “malla capilar” que cubre la totalidad del territorio urbano. Por el contrario, el Metro y BiciMAD muestran una participación más selectiva, concentrándose en corredores de alta capacidad y en la almendra central respectivamente, lo que refleja una jerarquía funcional clara entre modos de cobertura extensiva vs. intensiva.

7.3 Top 10 Zonas Más Versátiles

Interpretación del Ranking: Las zonas que encabezan este ranking exhiben valores de versatilidad cercanos a 1 (máxima entropía), lo que indica que no dependen de un solo modo de transporte. Estos nodos actúan como verdaderas “rótulas” del sistema, permitiendo a los usuarios cambiar de red (ej. de Metro a BiciMAD) con fluidez. Su importancia es sistémica: son los puntos donde se “cosen” las diferentes capas de la ciudad.

7.4 Análisis de Centralidad por Capa

Comparativa de Centralidad: Como se observa, la densidad y el grado medio varían drásticamente entre capas. La red de Autobuses presenta la mayor densidad de conexiones locales debido a la proximidad física de sus paradas, mientras que el Metro prioriza la eficiencia con menos conexiones pero más estratégicas. Esto valida numéricamente que cada capa tiene una topología optimizada para su función específica (cobertura local vs. velocidad a media distancia).

7.5 Robustez: Evolución de la Componente Gigante

Diagnóstico de Vulnerabilidad: Este gráfico es la evidencia más clara de la naturaleza “Scale-Free” del transporte madrileño. La divergencia extrema entre la línea de “Fallo Aleatorio” (que decae suavemente) y la de “Ataque Dirigido” (que colapsa la red eliminando apenas el 15% de los nodos) alerta sobre la fragilidad crítica de los hubs centrales. El sistema es robusto ante accidentes dispersos, pero necesita protección prioritaria en sus nodos vertebradores.

7.6 Distribución de Grados por Capa (Solo Nodos Activos)

Heterogeneidad de Grados: Los diagramas de violín revelan que, mientras la red de Autobuses tiene una distribución de grados más ancha y variable (reflejando zonas de muy alta y muy baja densidad), el Metro y el Parking muestran distribuciones más compactas. BiciMAD presenta un comportamiento bimodal interesante, sugiriendo una clara distinción entre estaciones periféricas del sistema y las estaciones centrales hiper-conectadas.

7.7 Relación Versatilidad vs Viajes

Desacoplamiento Estructural: La nube de puntos dispersa confirma la hipótesis de que infraestructura y demanda no siempre van de la mano. Existen zonas con alta versatilidad (eje Y alto) pero moderada demanda de viajes (eje X medio), que son infraestructuras estratégicas de conectividad más que destinos finales de masas. Esto sugiere una planificación que busca garantizar opciones de movilidad en toda la ciudad, no solo donde hay más viajeros.

7.8 Resumen Ejecutivo

Conclusiones Clave del Dashboard

Sistema Altamente Integrado: El 98.1% de las zonas ( 1528 de 1558 ) están conectadas en la componente gigante, demostrando una excelente cohesión del sistema.
Distribución por Modos: Bus ( 1346 zonas), Metro ( 228 ), BiciMAD ( 468 ), Parking ( 54 ), reflejando la jerarquía y capilaridad de cada modo.
Hubs Multimodales: Las 10 zonas con mayor versatilidad (≈1.0) actúan como intercambiadores críticos que integran múltiples modos de transporte.
Desacoplamiento Infraestructura-Demanda: La movilidad per cápita es uniforme (≈2 viajes/día) independientemente de la versatilidad, indicando planificación estratégica sobre demanda reactiva.
Vulnerabilidad Diferencial: Alta robustez ante fallos aleatorios pero fragilidad extrema ante ataques dirigidos a nodos centrales, especialmente en metro.
Topología Small-World: Clustering moderado (0.33) y alta conectividad global confirman una estructura óptima para transporte urbano.

7.8.1 Resumen Ejecutivo

El dashboard interactivo desarrollado permite explorar visualmente la estructura de la red multicapa de transporte de Madrid, revelando patrones estructurales clave del sistema urbano.

7.8.2 Hallazgos Principales

1. Alta Integración del Sistema El análisis revela que la red de Madrid está altamente integrada, con un 98.1% de las 1,558 zonas de intercambio conectadas en una única componente gigante. Esta elevada conectividad (solo 30 zonas aisladas) demuestra la efectividad de la planificación del transporte madrileño, garantizando que prácticamente toda la ciudad sea accesible mediante el sistema multimodal.

2. Identificación de Hubs Multimodales Críticos El análisis de versatilidad identifica 10 zonas estratégicas (IDs: 164, 446, 101, 41, 283, 56, 239, 814, 71, 534) con valores de entropía superiores a 0.99, indicando una presencia equilibrada de múltiples modos de transporte. Estas zonas actúan como verdaderos intercambiadores multimodales que vertebran el sistema.

3. Desacoplamiento entre Infraestructura y Demanda Individual Un hallazgo fundamental es que la versatilidad de las zonas (complejidad multimodal) no correlaciona con la intensidad de viajes per cápita, que se mantiene uniforme en aproximadamente 2 viajes/día en todos los distritos (rango: 1.99-2.13). Esto demuestra que la ubicación de los intercambiadores responde a criterios estratégicos de planificación de red y accesibilidad global, no a una mayor demanda local.

4. Topología de Red de Mundo Pequeño El coeficiente de clustering de 0.33 combinado con la alta conectividad global indica que Madrid presenta una topología de red de mundo pequeño (small-world network). Esta estructura, con alta conectividad local (clustering) y baja distancia global, es óptima para el transporte urbano, permitiendo desplazamientos eficientes con pocos transbordos.

5. Densidad Apropiada para Redes Espaciales La baja densidad observada (0.29%) es coherente y apropiada para redes de transporte urbano basadas en proximidad espacial. A diferencia de redes abstractas, las redes de transporte solo conectan nodos geográficamente próximos, resultando naturalmente en bajas densidades que no comprometen la funcionalidad del sistema.

7.8.3 Implicaciones para la Gestión del Sistema

Protección de Nodos Críticos: Las 10 zonas de mayor versatilidad identificadas deben considerarse críticas para la resiliencia del sistema y recibir atención prioritaria en mantenimiento y gestión de incidencias.
Planificación de Capacidad Uniforme: La homogeneidad en viajes per cápita sugiere que las inversiones en capacidad deben distribuirse equilibradamente entre distritos, sin concentrarse solo en intercambiadores.
Valor de la Multimodalidad: La alta versatilidad de ciertas zonas confirma el valor estratégico de invertir en intercambiadores multimodales que faciliten transbordos fluidos entre modos.

7.8.4 Nota sobre Metodología del Dashboard

El dashboard implementado utiliza visualizaciones interactivas (plotly, DT) que permiten explorar dinámicamente los datos. Las métricas de centralidad, versatilidad y robustez se calculan sobre la red agregada multicapa, proporcionando una visión holística del sistema que no sería visible analizando cada modo de transporte por separado.

La integración de datos de movilidad (viajes por distrito) con métricas estructurales de red permite conectar la topología de la infraestructura con los patrones reales de comportamiento de los usuarios, ofreciendo insights valiosos para la planificación urbana.

8 Conclusiones

En esta memoria se ha construido y analizado una red multicapas de transporte para Madrid utilizando exclusivamente los datos proporcionados:

Las paradas se han agrupado en zonas de intercambio mediante clustering espacial.
Se ha definido una red multicapa con capas para bus, metro, bicimad y parking.
Se ha calculado la centralidad de las zonas en cada capa y un indicador de versatilidad de grado que muestra qué zonas son importantes en varios modos a la vez.
Se ha generado una visualización multicapa con muxViz que hace explícita la estructura en capas.
Se ha estudiado la robustez estructural de la red agregada frente a fallos dirigidos en nodos centrales de la capa metro.

A partir de aquí, se podrían ampliar los análisis:

Incorporando información de líneas reales (no solo proximidad espacial).
Usando otras medidas de centralidad (betweenness, eigenvector, PageRank).
Implementando modelos de cascadas más cercanos a la literatura de redes interdependientes.

El código está preparado para ejecutarse directamente (asumiendo que los CSV y los paquetes necesarios están disponibles) y centrado únicamente en los datos reales de Madrid, tal y como se pedía en el enunciado.

8.1 Implementación Técnica y Uso de muxViz

Para la realización de este estudio, la librería muxViz en R ha sido la herramienta fundamental para el análisis matemático. El proceso se ha estructurado en fases clave:

Modelado de la Red: Hemos seguido la lógica matemática de muxViz para construir la Matriz Supra-Adyacente. Dado el tamaño de la red, hemos implementado esta construcción utilizando matrices dispersas optimizadas (Matrix package) en lugar de la función original BuildSupraAdjacencyMatrixFromExtendedEdgelist(), garantizando así la estabilidad y evitando desbordamientos de memoria, manteniendo intacta la definición del tensor supra-adyacente.
Análisis de Métricas: Para cuantificar la relevancia real de las paradas, hemos calculado:
- Grado Multicapa: Mediante GetMultiDegree(), obteniendo la conectividad total en el supra-grafo.
- PageRank Multicapa: Para este cálculo, hemos aplicado el algoritmo de PageRank directamente sobre el supra-grafo construido. Esta implementación es matemáticamente equivalente a GetMultiPageRankCentrality() pero se ejecuta mediante el motor optimizado en C++ de igraph, lo que asegura precisión numérica y velocidad de ejecución. Esta métrica es superior a la centralidad clásica, ya que identifica los nodos que son importantes por su capacidad de conectar flujos complejos a través de diferentes modos.
Visualización: Aunque muxViz proporciona funciones nativas para la visualización de redes multicapa, en este trabajo no se ha utilizado directamente la función original plot_multiplex3D(). Esto se debe a la existencia de errores de ejecución y problemas de compatibilidad con versiones recientes de R, documentados en los issues abiertos del repositorio oficial de muxViz en GitHub, que impedían su correcto funcionamiento.

Por este motivo, se ha empleado una versión modificada de dicha función, adaptada a partir del código original de muxViz, manteniendo su lógica de representación tridimensional y la filosofía del modelo multicapa. Esta adaptación permite generar una visualización 3D coherente en la que se aprecia claramente la estructura en capas (Bus, Metro, BiciMAD y Parking) y su alineación espacial, garantizando al mismo tiempo la estabilidad y reproducibilidad del análisis. Adicionalmente, la visualización se ha complementado con un dashboard interactivo en plotly para la exploración dinámica de las métricas multicapa.

Redes Multicapas en el Transporte de Madrid

Jairo Madrigal Cutillas, Mónica Sánchez Bellido, Nicolás Martín Morcillo, Eric Mas Belén

29/01/2026